Quantum spin chains in a magnetic field 
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We demonstrate that the "worm" algorithm allows very ef- 
fective and precise quantum Monte Carlo (QMC) simulations 
of spin systems in a magnetic field, and its auto-correlation 
time is rather insensitive to the value of H at low temperature. 
Magnetization curves for the s = 1/2 and s = 1 chains are 
presented and compared with existing Bethe ansatz and exact 
diagonalization results. From the Green function analysis we 
deduce the magnon spectra in the s = 1 system, and directly 
establish the "relativistic" form E(p) = (A 2 +v 2 p 2 ) 1/2 of the 
dispersion law. 



I. INTRODUCTION 

Quantum spin chains have recently been the subject of 
intensive theoretical, experimental and numerical stud- 
ies. Many methods were developed and applied first to 
these systems. In this paper we present precise calcula- 
tions of the magnetization curves, m z (H) = S z /L (where 
S z = Ylj—i s z(i) is the projection of the total spin on the 
direction of magnetic field, and L is the number of spins 
in the system), and excitation spectra in isotropic anti- 
ferromagnetic Heisenberg chains. 

The magnetization is probably the easiest quantity 
to measure experimentally, and, on the contrary, the 
most difficult to calculate numerically, as far as quantum 
Monte Carlo (QMC) methods are concerned. In canoni- 
cal ensemble (or ^-conserving) schemes, m z can be ob- 
tained only from the spin and field dependence of the 
ground state energies E(S Z ,H) witb_subsequent extrap- 
olation to the thermodynamic limit £1 Precise QMC cak 
culations of that kind require enormous numerical effort .a 
Also, most of the schemes rely-on the Suzuki- Trotter dis- 
cretization of imaginary time,B which introduces system- 
atic errors. On the other hand, the loop cluster update 
(LCU) algorithm is suffering from exponential slowing 
down in the most interesting cases (see Sec. IB due to 
small acceptance rates in a finite magnetic fieldD We are 
not aware of any large system simulations done using 
LCU in strong magnetic field. Usually precise calcula- 
tions of th&-.m z curves were done by exact diagonaliza- 
tion (ED)BEl and density matrix renormalization group 
(DMRG) methods^ 

With thfijdevelopment of the continuous-time "worm" 
algorithnxMa (it is called "worm" algorithm because the 



configuration space of the system is sampled through the 
worm- like motion of world line discontinuities, or, more 
generally, source operators) which works with the grand 
canonical ensemble (i.e., samples all S z ), the magnetiza- 
tion is calculated as easily as any other "standard" ther- 
modynamic quantity like the energy. Furthermore, in 
the same calculation one collects statistics for the Green 
function, G(r, r) = — T < s~(0, 0)s + (r, r) >, where T 
stands for the time ordering, and = s x ± is y . Thus 
in a single QMC simulation done at low temperature one 
can directly measure not just E(T, H), but also m z {H), 
the single-particle excitation spectrum, critical indices of 
G(i,r), the susceptibility x(-ff) = dm z {H)/dH, the spin 
stiffness A S (H), and more. This information may be used 
for precise calculations of spin gaps and magnon disper- 
sion curves. 

In what follows we first verify, in Sec. II, by perform- 
ing an autocorrelation times analysis, that the "worm" 
algorithm does not suffer strongly from slowing down in 
finite magnetic field and compare its performance with 
the continuous-time LCU method for s = 1/2 chains. In 
Sec. Ill we present our data for the magnetization curves 
in s = 1/2 and s = 1 chains and compare them with the 
results of ED and DMRG studies. In Sec. IV we calculate 
Haldane gaps and magnon dispersion curves for the s = 1 
system. From our data we confirm the existence of a 
shallow minimum in x(H) around H — 3.2 (we use units 
such that magnetic field, temperature, magnon velocity, 
and spin gaps are measured in terms of the Heisenberg 
exchange coupling constant J), found previously in the 
ED studies on small systems.^ We exclude the possibility 
that this minimum is a finite-size effect, or an artifact of 
extrapolation procedure used in Ref. |l|. Our result for 
the spin gap 



A(s = 1) = 0.4105(1) 



(1) 



shows the accuracy to which this quantity may be mea- 
sured in a single QMC simulation using the "worm" 
algorithm, and .-.agrees with the ED result A(s = 
1) = 0.41049(2)0 and the DMRG results A(s = 1) = 
0.41050(2), in Ref. | and A(s = 1) = 0.4104892(2) in 
Ref. [22[ The magnon dispersion curve fits perfectly to 
the relativistic form E 2 (p) — A 2 + v 2 p 2 , where v is the 
magnon velocity, and p is the magnon momentum. We 
deduced 



v(s = 1) = 2.48(1) (2) 

from the fit; this result is as accurate as the best known 
DMRG values deduced from the correlation length v(s = 
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FIG. 1. Illustration of the loop cluster update (LCU) algo- 
rithm on a spin dimer. 



FIG. 2. Illustration of the "worm" method on a spin dimer. 



1) = 2|475(5)EJ and magnetization studies v(s = 1) = 
2.49(1)H We also found, that for the s = 1 chain with 
L = 100 lattice sites even two-magnon states have vis- 
ible corrections to the hard-core bpaon picture which is 
sometimes used in fitting the data.lill 

II. CRITICAL SLOWING DOWN OF THE 
"WORM" AND LCU ALGORITHMS 

The LCU algorithm has been applied successfully to 
quantum spin systems near phase .transitions, without 
any sign of critical slowing down.E3~til As soon as a mag- 
netic field is turned on, however, the LCU updates show 
a severe, exponential slowing down. The reason is that in 
the standard LCU method the field cannot be taken into 
account when building the loop cluster, but is treated as 
a global weight, modifying the flipping probabilities of 
the loop clusters. 

We illustrate this problem in the simplest case, two 
spins in a magnetic field 

-ffdimcr = </SiS 2 — H(Si + S%). (3) 

At H — J the ground state of this dimer system is 
degenerate, with the S z = 1 triplet and the spin singlet 
both having energy —3/4 J. In the world line picture the 
S z = 1 triplet is represented by two straight world lines 
with Sf = +1/2 (Figs. |l|a and |a). Starting from this 
configuration the LCU method proposes to flip one of the 
world lines (Fig. gb). The new configuration however is 
energetically unfavorable as its energy is —1/4 J. It will 
be accepted only with a probability p = cxp(— /3J/2), 
which takes an exponential amount of time. Once this 
configuration has been accepted however it can quickly 
relax to the energetically favorable singlet state, in which 
a world line gains exchange energy by jumping between 
the two sites, as shown in Fig. [He. Returning from the 
singlet to the triplet state is again very unlikely, since the 
LCU first has to create a world line configurations with 
two straight world lines. Since there is a finite probabil- 
ity density (2/ J)dr for having the world line jump to the 
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FIG. 3. Integrated autocorrelation times for the magneti- 
zation on a s = 1/s dimer and a chain for both the loop 
cluster update (LCU) algorithm and the "worm" algorithm. 

neighboring site if the two spins are antiparallel, remov- 
ing all the jumps is again exponentially rare. Thus in the 
LCU for a spin- 1/2 dimer at the critical point J = H in 
a magnetic field the autocorrelation time for the magne- 
tization is exponentially large 

i3g£ocexp0H/2. (4) 

The dynamic exponent, defined as r cx f3 z is z u = oo. 

We measured the integrated autocorrelation time us- 
ing a multi cluster loop algorithm and confirmed the ex- 
ponential slowing down, as shown in Fig. |^. In higher 
dimensional systems the scaling remains exponential. 

The above discussion has pointed out that in going 
from the triplet to the singlet state we trade potential 
energy for exchange energy. The LCU method does this 
in two steps, first paying a huge loss in one, before gaining 
the other. The worm algorithm on the other hand can 
make these trades on a local basis, as it allows configu- 
rations with broken world lines, as is illustrated in Fig. 
pi. Starting from the triplet state the worm algorithm 
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flips only a short segment of the world line, with a high 
acceptance rate (Fig. ||b). Once this short segment has 
been created the world lines can gain exchange energy 
by jumping between the sites, as shown in Fig. |]c. The 
worm ends make a random walk along the time direction 
and need a time proportional to 1 to wind once around 
the time direction and thus change the magnetization. 
The autocorrelation time is 
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(5) 



and the dynamical exponent z™^™ = 1. The factor 2/3 
in the denominator normalizes the time by the minimum 
time needed in any algorithm, which is proportional to 
the space-time volume. 

We measured r in the worm algorithm and confirmed 
this linear scaling for the dimer. In the worm algorithm 
we measure the times in units of N (3 proposed local up- 
dates of the worm configuration, where N is the number 
of spins. As shown in Fig. ||, when (3H > 5 the worm 
algorithm is much faster than the LCU algorithm. 

In higher dimensions in an TV-spin system, the magne- 
tization fluctuates in general by \J N/ ' j3. To achieve an 
independent configuration one has to again spend time 



N/3 



(6) 



and thus z w 0. Close to a critical point the scaling 
is only slightly worse in one dimension. Modeling the 
magnons as hard-core bosons and assuming a dispersion 
e(k) oc k n the number of magnons at low temperatures 
is M (x Nf3 d / n , and the autocorrelation time 
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, 1) cx max^ 1 ^/™, const) (7) 



where d is the dimension. For a one-dimensional spin 
chain at either critical point, close to the spin gap H = 
A or close to the fully polarized state H = 2zJs the 
dispersion is quadratic (n — 2) and thus z — 1/2. In 
higher dimensions z w is predicted from this argument. 
Measurements of r on a chain, also shown in Fig. || 
confirm that r indeed increases only very slowly with 
0h. 

The above discussion assumed that the source opera- 
tors in the worm algorithm perform a random walk on 
the lattice. In the simple worm algorithm however they 
do not perform a pure random walk, but are guided by 
the Green function. The worm movements can be bi- 
ased by the inverse of the Green function (or an estimate 
thereof) though, which makes them a random walk. 

Finally we wish to note that work is in progress on 
improvingj-the scaling of the loop algorithm in a mag- 
netic field.ll3 For general couplings in the dimer case and 
for ferromagnetic couplings on any lattice the exponen- 
tial slowing down can be removed. However in the most 
interesting case, discussed here, an antiferromagnet on 
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FIG. 4. Magnetization curve for the s = 1/2 chain. The 
QMC data are shown by filled circles and the Bethe ansatz 
data by open circles (the point sizes is exaggerated, especially 
for the Bethe ansatz data, which were calculated up to four 
meaningful digits. 



a chain or higher dimensional lattice, the slowing down 
remains exponential, although the prefactor can be re- 
duced. 



III. MAGNETIZATION CURVES 

We have seen in the previous Sec. II that the "worm" 
algorithm is essentially insensitive to the value of H . 
As a testing case we computed the magnetization curve 
m z (H) for the s = 1/2 chain with L — 100 spins and 
verified that it agrees with the available Bethe ansatz 
solution.113 Within our error bars the QMC and Bethe 
ansatz data arc in general indistinguishable when calcu- 
lated at the same value of magnetic field. At (3 = 80 and 
very small fields H one can however see in Fig. ^ the typi- 
cal finite-size oscillations, which appear when (3 ~ 2nv/ L 
and only few magnons are present in the system; at T = 
these oscillations convert into a stair-case curve. These 
effects can be reduced by going to larger L. The whole 
curve takes about 40 hours of CPU time on an HP-UX 
9000/735. 

The case of s — 1 is more intriguing. A Bethe ansatz 
solution is not available in this case, and, as mentioned 
above, one has to rely on ED or DMRG. The magneti- 
zation curve for the s = 1 chain was calculated in Ref. 
[IJ. Some aspects of this study may rise questions. Since 
the largest system size was only L = 16, the expected 
square-root singularity near H c \ = A(s = 1) was hardly 
seen, although it is possible do recoverpthis singularity 
using an appropriate finite-size analysis. 113 

Also, after extrapolating results for E(S Z ,H) to the 
thermodynamic limit to deduce m z {H) and subsequent 
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FIG. 5. Magnetization curve of the s = 1 antiferromagnetic 
Heisenberg chain. 

differentiation of m z (H), a very shallow minimum in 
x(H) around H ss 3.2 was predicted. The amplitude 
of the x(-ff) variation around the minimum was only a 
few percent, and one might suspect that the procedure 
of eliminating finite-size effects was simply not accurate 
enough. 

In Fig. ^|we show our data for m z (H) calculated with 
high accuracy to ensure reliable values for the field deriva- 
tive x(H)- Obviously, there is a wide region in H where 
m z {H) is almost perfectly linear. In Fig. || we demon- 
strate that, beyond any doubts, there is a minimum in 
x(-ff) around H « 3.2, and that the technique developed 
in Ref. ^ for eliminating finite-size effects in small clus- 
ters is remarkably accurate once the system size is larger 
than ilie correlation length (£ = 6.03(1) for the s = 1 
chainB). The origin of this minimum is not understood 
theoretically yet. Furthermore, we calculated x(-ff) in 
the region between H = 2.8 and H = 3.2, where dx/dH 
is negative, at higher temperature. If the minimum was 
any kind of finite size effect due to level quantization, it 
would have been strongly affected by temperature varia- 
tion. We see, that higher-temperature points, marked by 
crosses in Fig. are not affected within error bars, and 
also give a negative dx/dH in this region. 

Clearly, for the system with L = 100, we already see 
the square-root singularity m z ~ y/H — A near the spin 
gap values. It is smeared out very close to the gap due to 
finite temperature. The physics of m z (H) behavior near 
the gap is usually described by a free-fermion model with 
periodic/antiperiodic bopadary conditions depending ©a. 
the odd/even value of S z ^3 or a hard-core boson model,t2l 



FIG. 6. Derivative of the magnetization with respect to the 
field x(H), for the s = 1 antiferromagnetic Heisenberg chain 
at P = 100/ J. The region around the minimum is shown 
enlarged in the inset. Higher temperature data (at j3 = 50/ J) 
are marked by crosses. 




FIG. 7. Magnetization curve of the s = 1 antiferromagnetic 
Heisenberg chain close to the gap A on an L = 100 site chain 
at P = 100. The solid curve shows the prediction of the 
hard-core boson theory assuming A = 0.4105 and v = 2.48. 
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which is exact in the dilute limit m z — > 0. It is tempting 
then to fit the QMC data near H c \ to the prediction of 
the hard-core boson theory at finite, but very low tem- 
perature. In Fig. we show the result of such a fit, done 
using the following ansatz for the magnon spectrum 

E(p) = VA 2 + v V , (8) 

with A(s = 1) = 0.4105 and v(s = 1) = 2.48. Since finite 
size corrections in a chain with periodic boundary condi- 
tions are exponentially small, we used the known values 
of A(s = 1) and v(s = 1). Surprisingly, the fit is not good 
already for m z > 0.02, i.e., when only two magnons are 
inserted into the system. Also, the theory predicts visible 
finite-size oscillations in x(-ff) at (3 = 100, which are not 
present in the calculated m z (H). Although it is possible 
to get an almost perfect fit of the QMC data to the hard- 
core boson theory up to m ~ 0.06 using v(s = 1) = 2.25, 
we rather suggest that the difference between the theo- 
retical and calculated m z is due to magnon-magnon in- 
teractions other than on-site repulsion. The gap and the 
magnon velocity were deduced independently with high 
accuracy from the Green function (see next Section) and 
exclude the possibility that v = 2.25. It means that one 
must be extremely cautious in extracting spin velocity 
from multi-magnon states, i.e. states with S z > 3, since 
longer range interactions between the magnons are usu- 
ally ignored in such studies lid 

We note, that the magnetization curve itself in Fig. 
[j] is sufficient to deduce the gap value rather accurately 
A(s = 1) = 0.410(2) even without the Green function 
analysis to which we proceed now. 

IV. MAGNON SPECTRA 

Since the "worm" algorithm collects "by passing" the 
complete histogram for the system's Green function, one 
may use it in a very efficient way to deduce the single- 
particle (magnon) spectrum. In the present study we 
demonstrate how well the method works for the case 
when an analytic continuation of G(i, r) to real frequen- 
cies is not necessary. 

At very low temperature and \H\ < H c i, the statistics 
is dominated by the ground-state configurations, and a 
piece of an extra world line, corresponding to the S z — 
±1 states. If H c i-H < H c i and T < H c i~H, then most 
of the G{p, t) histogram (in momentum representation) 
is determined by the S z = +1 configurations except for 
very short ~ 1 negative times, (since the typical length 
in time of an extra worldline is of order 1/(H C \ — H), 
we keep H rather close to the critical field in order to 
make the worldline longer; this allows one more precise 
determination of the Green function decay in time). In 
our simulations we used (3 — 200 and H c i — H ~ 10//3. 
One then expects that (counting E(p) from H) 

G(p,r)~e- E ^, (9) 
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FIG. 8. The magnon spectrum in the s = 1 chain with 
L — 100 spins calculated at j3 = 200/ J. The solid line is the 
fit to the "relativistic" spectrum, and the dotted line is the 
parabolic approximation with the same parameters for the 
gap and velocity. The zero-momentum result is also shown in 
the insert. 

i.e., a simple exponential decay from which the magnon 
spectrum is readily obtained as 

= HG(p,t)/G{ P ,t q )} ^ 

T - T 

This definition is normalization free, although normal- 
ization is not a problem and is fixed by the condition 
-G(i = 0,T = +0) = s + l. 

First of all, we verify that the decay of G(p, r) in time 
is purely exponential for r» 1 until the value of G be- 
comes too small and statistical fluctuations start domi- 
nating. We then use Eq. ([!(]) to determine the spectrum. 
For each momentum state the reference point was set to 
t = l/E(p), and then the value of E(p) was obtained 
from the best exponential fit to the data between To and 
Tmax, where T max (j>) was close to the onset of statistical 
fluctuations in G(p, r) (typically E(p)T max ~ 4 — 5). 
The error bars are estimated from the fluctuations of re- 
sults of up to twelve independent simulations. Our re- 
sults for the chain s = 1 are shown in Fig. ||. For the 
s = 1 chain we studied a system with L = 100 spins. The 
precision of the data, obtained in only 50 hours of CPU 
time on an HP-UX 9000/735, is sufficient to determine 
A(s = 1) = 0.4105(1) — the most accurate QMC result 
so far. No finite size corrections are necessary since these 
are exponentially small, given the rather short correla=j 
tion length jfht-s = 1. This value agrees with the EDeI 
and DMRG&E3 data. The magnon velocity was found 
to equal v(s = 1) = 2.48(1). This is also the most ac- 
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curate QMC result and a gree s with the DMRG values 
v(s = 1) = 2.475(5) in Ref. [l(], and v(s = 1) = 2.49(1) in 
Ref. 0. 

One note is in order here. From Fig. || we see, that 
the parabolic expansion E(p) — A + v 2 p 2 /2A is not ac- 
curate for the chain L = 100, except maybe for the first 
nonzero-momentum state, but the data are consistent 
with the form (||) for all the reliably determined states. 
This unambiguously confirms the "relativistic" ansatz for 
the dispersion law. The spin velocity is sometimes deter- 
mined from the ground-state energies E(S Z , L) for chains 
of comparable sizes assuming parabolic expansion and ig- 
noring magHon-magnon interactions (other than on-site 
repulsion) Jiil Since we found that both aspects have no- 
ticeable corrections to the data in chains with L ~ 100, 
we suspect that these corrections somehow compensate 
each other and the net result turns out to come right, 
e.g., the DMRG result obtained in Ref. [tl] is v = 2.49(1). 
Our data give the dispersion law directly and do not 
rely on any assumptions about its form. Also, they are 
not affected by the magnon-magnon interactions since 
the parameters H and T are chosen in such a way that 
only single-magnon virtual states are contributing to the 
statistics. 



V. CONCLUDING REMARKS 

In the present study we have shown that "worm" al- 
gorithm is very efficient for simulations of spin systems 
in finite magnetic field. Its autocorrelation time is al- 
most independent of (3 and H. For the isotropic anti- 
ferromagnetic Heisenberg chains s = 1/2 and s = 1 the 
accuracy of our data for the magnetization curves, spin 
gaps and magnon velocities is comparable with the Bethe 
ansatz, exact diagonalization and DMRG results. More- 
over, from these data we found visible longe-range (apart 
from onsite hard-core repulsion) magnon-magnon inter- 
actions, and unambiguously verify the relativistic form of 
the magnon dispersion law. Our results for the magnon 
spectrum are the most straightforward ones, i.e., they 
are not affected by the finite-size effects (L/£ > 10), spin 
excitations at the chain ends, or magnon-magnon inter- 
actions. Also, the numerical effort is not extreme, since 
one has to perform only one simulation to get the above 
results. This is made possible through the evaluation of 
the Green function of the system. 

It was demonstrated recently how to implement 
the "woxm" idea within the framework of the LCU 
method.EEl. Since A(s = 2) is rather small ~ 0.09 (see 
Refs. [2l|-p3|), one may deduce the gap and the spec- 
trum very accurately from the Green function even at 
H = making .use of improved estimators available in 
LCU algorithms^. 
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